home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 26 / Cream of the Crop 26.iso / os2 / octa209s.zip / octave-2.09 / libcruft / lapack / zunglq.f < prev    next >
Text File  |  1996-07-19  |  6KB  |  208 lines
  1.       SUBROUTINE ZUNGLQ( M, N, K, A, LDA, TAU, WORK, LWORK, INFO )
  2. *
  3. *  -- LAPACK routine (version 2.0) --
  4. *     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
  5. *     Courant Institute, Argonne National Lab, and Rice University
  6. *     September 30, 1994
  7. *
  8. *     .. Scalar Arguments ..
  9.       INTEGER            INFO, K, LDA, LWORK, M, N
  10. *     ..
  11. *     .. Array Arguments ..
  12.       COMPLEX*16         A( LDA, * ), TAU( * ), WORK( LWORK )
  13. *     ..
  14. *
  15. *  Purpose
  16. *  =======
  17. *
  18. *  ZUNGLQ generates an M-by-N complex matrix Q with orthonormal rows,
  19. *  which is defined as the first M rows of a product of K elementary
  20. *  reflectors of order N
  21. *
  22. *        Q  =  H(k)' . . . H(2)' H(1)'
  23. *
  24. *  as returned by ZGELQF.
  25. *
  26. *  Arguments
  27. *  =========
  28. *
  29. *  M       (input) INTEGER
  30. *          The number of rows of the matrix Q. M >= 0.
  31. *
  32. *  N       (input) INTEGER
  33. *          The number of columns of the matrix Q. N >= M.
  34. *
  35. *  K       (input) INTEGER
  36. *          The number of elementary reflectors whose product defines the
  37. *          matrix Q. M >= K >= 0.
  38. *
  39. *  A       (input/output) COMPLEX*16 array, dimension (LDA,N)
  40. *          On entry, the i-th row must contain the vector which defines
  41. *          the elementary reflector H(i), for i = 1,2,...,k, as returned
  42. *          by ZGELQF in the first k rows of its array argument A.
  43. *          On exit, the M-by-N matrix Q.
  44. *
  45. *  LDA     (input) INTEGER
  46. *          The first dimension of the array A. LDA >= max(1,M).
  47. *
  48. *  TAU     (input) COMPLEX*16 array, dimension (K)
  49. *          TAU(i) must contain the scalar factor of the elementary
  50. *          reflector H(i), as returned by ZGELQF.
  51. *
  52. *  WORK    (workspace/output) COMPLEX*16 array, dimension (LWORK)
  53. *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
  54. *
  55. *  LWORK   (input) INTEGER
  56. *          The dimension of the array WORK. LWORK >= max(1,M).
  57. *          For optimum performance LWORK >= M*NB, where NB is
  58. *          the optimal blocksize.
  59. *
  60. *  INFO    (output) INTEGER
  61. *          = 0:  successful exit;
  62. *          < 0:  if INFO = -i, the i-th argument has an illegal value
  63. *
  64. *  =====================================================================
  65. *
  66. *     .. Parameters ..
  67.       COMPLEX*16         ZERO
  68.       PARAMETER          ( ZERO = ( 0.0D+0, 0.0D+0 ) )
  69. *     ..
  70. *     .. Local Scalars ..
  71.       INTEGER            I, IB, IINFO, IWS, J, KI, KK, L, LDWORK, NB,
  72.      $                   NBMIN, NX
  73. *     ..
  74. *     .. External Subroutines ..
  75.       EXTERNAL           XERBLA, ZLARFB, ZLARFT, ZUNGL2
  76. *     ..
  77. *     .. Intrinsic Functions ..
  78.       INTRINSIC          MAX, MIN
  79. *     ..
  80. *     .. External Functions ..
  81.       INTEGER            ILAENV
  82.       EXTERNAL           ILAENV
  83. *     ..
  84. *     .. Executable Statements ..
  85. *
  86. *     Test the input arguments
  87. *
  88.       INFO = 0
  89.       IF( M.LT.0 ) THEN
  90.          INFO = -1
  91.       ELSE IF( N.LT.M ) THEN
  92.          INFO = -2
  93.       ELSE IF( K.LT.0 .OR. K.GT.M ) THEN
  94.          INFO = -3
  95.       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
  96.          INFO = -5
  97.       ELSE IF( LWORK.LT.MAX( 1, M ) ) THEN
  98.          INFO = -8
  99.       END IF
  100.       IF( INFO.NE.0 ) THEN
  101.          CALL XERBLA( 'ZUNGLQ', -INFO )
  102.          RETURN
  103.       END IF
  104. *
  105. *     Quick return if possible
  106. *
  107.       IF( M.LE.0 ) THEN
  108.          WORK( 1 ) = 1
  109.          RETURN
  110.       END IF
  111. *
  112. *     Determine the block size.
  113. *
  114.       NB = ILAENV( 1, 'ZUNGLQ', ' ', M, N, K, -1 )
  115.       NBMIN = 2
  116.       NX = 0
  117.       IWS = M
  118.       IF( NB.GT.1 .AND. NB.LT.K ) THEN
  119. *
  120. *        Determine when to cross over from blocked to unblocked code.
  121. *
  122.          NX = MAX( 0, ILAENV( 3, 'ZUNGLQ', ' ', M, N, K, -1 ) )
  123.          IF( NX.LT.K ) THEN
  124. *
  125. *           Determine if workspace is large enough for blocked code.
  126. *
  127.             LDWORK = M
  128.             IWS = LDWORK*NB
  129.             IF( LWORK.LT.IWS ) THEN
  130. *
  131. *              Not enough workspace to use optimal NB:  reduce NB and
  132. *              determine the minimum value of NB.
  133. *
  134.                NB = LWORK / LDWORK
  135.                NBMIN = MAX( 2, ILAENV( 2, 'ZUNGLQ', ' ', M, N, K, -1 ) )
  136.             END IF
  137.          END IF
  138.       END IF
  139. *
  140.       IF( NB.GE.NBMIN .AND. NB.LT.K .AND. NX.LT.K ) THEN
  141. *
  142. *        Use blocked code after the last block.
  143. *        The first kk rows are handled by the block method.
  144. *
  145.          KI = ( ( K-NX-1 ) / NB )*NB
  146.          KK = MIN( K, KI+NB )
  147. *
  148. *        Set A(kk+1:m,1:kk) to zero.
  149. *
  150.          DO 20 J = 1, KK
  151.             DO 10 I = KK + 1, M
  152.                A( I, J ) = ZERO
  153.    10       CONTINUE
  154.    20    CONTINUE
  155.       ELSE
  156.          KK = 0
  157.       END IF
  158. *
  159. *     Use unblocked code for the last or only block.
  160. *
  161.       IF( KK.LT.M )
  162.      $   CALL ZUNGL2( M-KK, N-KK, K-KK, A( KK+1, KK+1 ), LDA,
  163.      $                TAU( KK+1 ), WORK, IINFO )
  164. *
  165.       IF( KK.GT.0 ) THEN
  166. *
  167. *        Use blocked code
  168. *
  169.          DO 50 I = KI + 1, 1, -NB
  170.             IB = MIN( NB, K-I+1 )
  171.             IF( I+IB.LE.M ) THEN
  172. *
  173. *              Form the triangular factor of the block reflector
  174. *              H = H(i) H(i+1) . . . H(i+ib-1)
  175. *
  176.                CALL ZLARFT( 'Forward', 'Rowwise', N-I+1, IB, A( I, I ),
  177.      $                      LDA, TAU( I ), WORK, LDWORK )
  178. *
  179. *              Apply H' to A(i+ib:m,i:n) from the right
  180. *
  181.                CALL ZLARFB( 'Right', 'Conjugate transpose', 'Forward',
  182.      $                      'Rowwise', M-I-IB+1, N-I+1, IB, A( I, I ),
  183.      $                      LDA, WORK, LDWORK, A( I+IB, I ), LDA,
  184.      $                      WORK( IB+1 ), LDWORK )
  185.             END IF
  186. *
  187. *           Apply H' to columns i:n of current block
  188. *
  189.             CALL ZUNGL2( IB, N-I+1, IB, A( I, I ), LDA, TAU( I ), WORK,
  190.      $                   IINFO )
  191. *
  192. *           Set columns 1:i-1 of current block to zero
  193. *
  194.             DO 40 J = 1, I - 1
  195.                DO 30 L = I, I + IB - 1
  196.                   A( L, J ) = ZERO
  197.    30          CONTINUE
  198.    40       CONTINUE
  199.    50    CONTINUE
  200.       END IF
  201. *
  202.       WORK( 1 ) = IWS
  203.       RETURN
  204. *
  205. *     End of ZUNGLQ
  206. *
  207.       END
  208.